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We derive and investigate the microscopic model of the quantum magnet BiCu2P06 using band 
structure calculations, magnetic susceptibility and high-field magnetization measurements, as well 
as Exact Diagonalization (ED) and Density-Matrix Renormalization Group (DMRG) techniques. 
The resulting quasi-one-dimensional spin model is a two-leg antiferromagnetic ladder with frustrat- 
ing next-nearest-neighbor couplings along the legs. The individual couplings are estimated from 
band structure calculations and by fitting the magnetic susceptibility with theoretical predictions, 
obtained using full diagonalizations. The nearest-neighbor leg coupling Ji, the rung coupling J4, 
and one of the next-nearest-neighbor couplings J2 amount to 120 — 150 K, while the second next- 
nearest-neighbor coupling is J2 ~ J2/2. The spin ladders do not match the structural chains, and 
although the next-nearest-neighbor interactions J2 and J2 have very similar superexchange path- 
ways, they differ substantially in magnitude due to a tiny difference in the 0-0 distances and in the 
arrangement of non-magnetic PO4 tetrahedra. An extensive ED study of the proposed model pro- 
vides the low-energy excitation spectrum and shows that the system is in the strong rung coupling 
regime. The strong frustration by the next-nearest-neighbor couplings leads to a triplon branch 
with an incommensurate minimum. This is further corroborated by a strong-coupling expansion up 
to second order in the inter-rung coupling. Based on high-field magnetization measurements, we 
estimate the spin gap of A ~ 32 K and suggest the likely presence of antisymmetric Dzyaloshinskii- 
Moriya anisotropy and inter-ladder coupling J3. We also provide a tentative description of the 
physics of BiCu2P06 in magnetic field, in the light of the low-energy excitation spectra and numer- 
ical calculations based on ED and DMRG. In particular, we raise the possibility for a rich interplay 
between one- and two-component Luttinger liquid phases and a magnetization plateau at 1/2 of the 
saturation value. 



PACS numbers: 75.50.-y, 75.30.Et, 75.10.Jm, 71.20.Ps 

I. INTRODUCTION 

One-dimensional (ID) spin systems are in the fo- 
cus of the present-day research due to a range of un- 
usual low-temperature properties governed by quantum 
effects. The primary ID spin model is the uniform spin-i 
Heisenberg chain that has a peculiar gapless excitation 
spectrumii Numerous model compounds and the large 
set of theoretical tools in one dimension made extensive 
comparisons between experiment and theory possible: for 
example, the universal scaling of spin excitations in the 
uniform spin-i Heisenberg chain was proposed theoreti- 
cally and later confirmed experimentally. 2 A number of 
studies successfully extended the model by including in- 
terchain couplings and discussed the trends for the or- 
dering temperature depending on the topology and mag- 
nitude of interchain couplings,— - — 

Alterations in the chain topology lead to a dramatic 
change in the magnetic properties. For example, there 
are several options to switch from the gapless spectrum 
of the uniform spin-^ chain to a gapped spectrum. The 
latter offers an exciting opportunity to close the spin gap 
by an external magnetic field and to observe unusual 
phenomena, such as Luttinger liquid (LL) physics and 
the Bose-Einstein condensation of triplons in the gap- 



less high-field phased The simplest way to introduce a 
spin gap into a ID system is to alternate the exchange 
couplings along the chain.— Another option is the frus- 
tration of the chain by next-nearest-neighbor couplings. 8 
Finally, several chains can be joined into a spin ladder 
that shows a spin gap for an even number of legs£ De- 
spite the relatively simple chain geometries, such models 
are rather difficult to realize experimentally. There is 
still no experimental observation of the LL phase in the 
alternating spin-i chain, and experimental examples of 
gapped frustrated spin chains are rare^ The quest for 
spin-ladder systems was more successful. For example, 
recently a remarkable mapping of high-field properties 
onto the LL model in a (CsH^N^CuB^ compound was 
performed j 1 - 1 - - — 

Combining different features of the modified chain 
topology (alternation, frustration, and coupling into 
a ladder), one can achieve further interesting proper- 
ties. For example, frustrated spin chains with alternat- 
ing nearest-neighbor couplings are predicted to exhibit 
a magnetization plateau for a certain range of model 
parameters* 1 ^ However, this prediction has never been 
tested experimentally due to the lack of proper model 
compounds. The problems with finding experimental 
realizations of certain spin models call for an alterna- 
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tive approach: the investigation of complex ID models, 
stimulated by real materials. In the following we show 
that the recently discovered spin-^ compound BiCu2POg 
closely corresponds to an interesting quasi-lD spin model 
combining all the three aforementioned features: frus- 
tration, spin-ladder geometry, and alternation of next- 
nearest-neighbor exchange couplings. 

Despite previous experimental and computational 
studies^ - — the microscopic model of BiCu 2 PC>6 remains 
controversial. To resolve this controversy, we apply a 
range of state-of-the-art computational techniques that 
reveal an accurate spin model and allow for a precise com- 
parison with the experimental results. First, we analyze 
the crystal structure and outline the previous reports in 
Sec. [11] After a brief description of the methods fSec. lIII[) . 
we proceed to extensive band structure calculations, de- 
rive a consistent spin model, and discuss the non-trivial 
implementation of this model in the crystal structure of 
BiCu2PC>6 fSec. IIV|) . In Sec. [V] we report the magnetic 
susceptibility and the high-field magnetization measure- 
ments that challenge the proposed spin model and un- 
ambiguously measure the spin gap. Finally, we perform 
model simulations, investigate the microscopic physics of 
BiCu2PC>6 at low energies (Sec. |VT|) and in the presence 
of magnetic field fSec. lVII|) . and conclude our study with 
a brief discussion and summary in Sec. IVIII1 



II. CRYSTAL STRUCTURE AND MAGNETIC 
PROPERTIES 

The crystal structure of BiCu2P06 (Fig. [I]) shows pro- 
nounced ID features with complex ribbons running along 
the b direction. 18 Each ribbon is formed by dimers of 
edge-sharing CUO4 plaquettes. The plaquettes of the 
neighboring dimers share corners (oxygen sites), while 
the ncxt-nearcst-neighbor dimers are additionally con- 
nected by PO4 tetrahedra. The spatial arrangement of 
the magnetic Cu atoms features both the spin-ladder and 
frustrated-spin-chain geometries (see Figs. Q]and|4]). The 
stacking of the dimers reminds of the spin ladder with 
the leg coupling J\ and the rung coupling Yet, 
the interactions J\ follow a zigzag pattern and form a 
frustrated spin chain, once the couplings between next- 
nearest neighbors are considered. The situation is fur- 
ther complicated by the two inequivalent Cu positions, 
leading to inequivalent next-nearest-neighbor couplings 
J 2 andJ 2 3 

The complex crystal structure of BiCu2PC>6 led to a 
controversy regarding the appropriate spin model of this 
compound. Koteswararao et alr^ emphasized the spin- 
ladder feature of the structural ribbons and considered 
BiCu2PC"6 as a system of J\ — J3 ladders that are cou- 
pled by the inter-ribbon interaction J4. This interpreta- 
tion prevailed in further studies, focused on the effects 
of dopingi 17 ' 21 ~ — However, band structure calculations, 
reported by the same authors?^ clearly showed sizable 
next-nearest-neighbor couplings J2 and J' 2 that would in- 




FIG. 1. (Color online) Crystal structure of BiCu2P06 with 
ribbons comprising CUO4 plaquettes and PO4 tetrahedra 
(top) and the spin model (bottom). Open and shaded cir- 
cles denote the two inequivalent Cu positions, while the larger 
dark circles label the Bi atoms. More details on the structure 
are shown in Fig. [J] The model is of the spin-ladder type 
and comprises four inequivalent couplings: the leg coupling 
Ji, the rung coupling J4, and the frustrating next-nearest- 
neighbor leg couplings J2 and J'i- Note that the two legs of 
the ladder reside on different structural ribbons. 



evitably frustrate the system. 

Although similar at a first glance, Mentre et al*& sug- 
gested a somewhat different spin model. Using inelas- 
tic neutron scattering (INS) and band structure calcu- 
lations, they showed that the ladders are formed by the 
couplings Ji and J4, while the intra-ribbon interaction J3 
is an inter-ladder coupling. To fit the INS data, Mentre 
et al. also had to include the next-nearest-neighbor cou- 
pling J 2 , but the difference between J 2 and J 2 could not 
be resolved. 

Experimentally, BiCu2POg is a spin-gap material with 
a singlet ground state (no long-range ordering). The 
substitution of Cu by non-magnetic Zn atoms destroys 
the spin gap and leads to a spin freezing ! 17 i 21 These fea- 
tures are fairly general and can be assigned to a range 
of simple ID spin models (alternating chain, frustrated 
chain, two- leg ladder). However, the experimental data 
can not be described well by any of these models (see also 
Sec. |Vj. The previous report a 15 ' 16 evidence the combi- 
nation of the ladder-type geometry and the frustration 
by next-nearest-neighbor couplings. Yet, the precise way 
of this combination and, more importantly, the resulting 
physics remain unclear. 



III. METHODS 

To evaluate the individual exchange couplings in 
BiCu2POg, we performed scalar-relativistic density func- 
tional theory (DFT) band structure calculations using 
the full-potential local-orbital FPLO code (version 8.00- 
31) ,24 The calculations were done in the framework of 
the local (spin) density approximation [L(S)DA], employ- 
ing the exchange-correlation potential by Perdew and 
Wangj2£ The symmetry-irreducible part of the first Bril- 
louin zone was sampled by a mesh of 512 fc-points for the 
crystallographic unit cell and 64 k points for the super- 
cells. 

Superexchange couplings in insulating Cu +2 com- 
pounds are intimately related to strong electronic corre- 
lations that cannot be properly treated within L(S)DA. 
To account for the correlation effects, we used two ap- 
proaches. First, we mapped the half-filled LDA Cu 3d 
bands via an effective one-band tight-binding (TB) model 
onto a Hubbard model. Then, antiferromagnetic (AFM) 
exchange integrals were derived from the expression of 
the second-order perturbation theory. This procedure 
is referred below as the model approach. In the second 
(supercell) approach, the correlation effects were treated 
in a mean-field approximation within the band structure 
calculations by applying the LSDA+C/ method. 26 The 
on-site Coulomb repulsion parameter Ud was varied in 
the 6 — 8 eV range while the on-site exchange pa- 
rameter Jd was fixed to 1 eV. Total energies for different 
types of collinear magnetic ordering were obtained within 
the crystallographic unit cell and the two supercells, dou- 
bled along the b or c directions. The calculated energies 
were mapped onto a Heisenberg model, and individual 
exchange couplings were derived. More details on the 
computational procedure are given in Sec. IIVI 

The resulting spin model was compared to the exper- 
imental results from magnetic susceptibility and high- 
field magnetization measurements. Powder samples of 
BiCu2POg were prepared by firing a stoichiometric mix- 
ture of Bi 2 3 (99.9 % purity), CuO (99.99 % purity), 
and NH4H2PO4 (99.9 % purity) in air. The mixtures 
were first annealed at 400 °C for 10 hours and then at 
850 °C for 40 hours with one intermediate grinding. The 
resulting samples were single-phase, as confirmed by x- 
ray diffraction (STOE STADI-P diffractometer, CuK Qi 
radiation, transmission geometry) . The magnetic suscep- 
tibility was measured in fields up to 5 T in the temper- 
ature range 2 — 700 K using a Quantum Design MPMS 
SQUID magnetometer. 

High-field magnetization measurements were per- 
formed at Hochfeld-Magnetlabor Dresden at 1.4 K tem- 
perature in fields up to 60 T using a pulsed magnet. De- 
tails of the measurement technique are given in Ref. l3lL 
The curves measured on increasing and decreasing field 
coincided, indicating the lack of any irreversible effects 
upon magnetization of the sample. 

Thermodynamic properties of the BiC^POe spin 
model were calculated by a full diagonalization for fi- 




FIG. 2. (Color online) Total and site-projected DOS obtained 
from LDA. The vertical line at zero energy denotes the Fermi 
level Ef- The bands near the Fermi level primarily comprise 
Cu and O states. The shading in the plot denotes the Cu-3d 
states, while the dashed line represents the O 2p states. 

nite lattices with N = 16 and 20 sites and periodic 
boundary conditions. To obtain the low-energy excita- 
tions, we performed Exact Diagonalizations (ED) using 
the Lanczos algorithm that allowed to extend the sys- 
tem size up to N = 36. The results are well converged 
with respect to the system size even for N = 16 and 
20, thus the finite-size effects for the spin model un- 
der consideration are relatively small. To obtain the 
magnetization process of BiCu 2 POg, we have used, in 
addition to ED, the Density Matrix Renormalization 
Group (DMRG) 32 ' 33 method with open boundary con- 
ditions with up to 128 rungs. Further details are given 
in Sees. El ED andEED 

IV. DERIVATION OF THE SPIN MODEL 

Spin models with exchange couplings derived from 
DFT have been previously reported in Refs. [l5| and [l6l . 
However, the analysis remains incomplete, since the two 
inequivalent next-nearest-neighbor couplings (between 
crystallographically different Cu sites) were considered to 
be equivalent. In the following, we apply two complemen- 
tary approaches that evaluate all the relevant exchange 
integrals and establish the microscopic model. Addition- 
ally, we analyze in detail the structural features that 
cause the unusual implementation of the ladder-type spin 
lattice in BiCu2P06- 

A. LDA and model approach 

Fig. [2] shows the LDA density of states (DOS) of 
BiCu2POg. The valence band spectrum is formed mainly 
by copper 3d and oxygen 2p orbitals, with a sizable con- 
tribution from phosphorous 3p orbitals below —3 eV. 
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FIG. 3. (Color online) LDA band structure (thin blue 
lines), the WF-based fit of the tight-binding model (bright 
orange dots), and the contribution of the Cu d x 2_ y 2 or- 
bital (dark purple dots). The high-symmetry fc-path 
in terms of the reciprocal lattice parameters is as fol- 
lows: r(0,0,0), X(0.5,0,0), 5(0.5,0.5,0), F(0, 0.5,0), L, 
2(0,0,0.5), [7(0.5,0,0.5), #(0.5,0.5,0.5), T(0, 0.5, 0.5). The 
bands are highly dispersive along X — S, Y — L — Z , and U — R 
which represent the leading interactions within the crystallo- 
graphic be plane and the quasi-2D nature of the system. 



The states above —0.6 eV are formed by the Cu 3d x a_ y a 
orbital, in agreement with the expected ligand-ficld 
splitting^ The shapes and positions of the bands close 
to the Fermi level (Ep) are somewhat different from the 
iV th -order muffin-tin orbital (NMTO) result of Ref. GJ, 
where the Cu 3d x 2_ y 2 bands are separated from the 
lower- lying bands. This difference represents a known 
shortcoming of the NMTO method^ To check our find- 
ings, we repeated the calculation using the full-potential 
code Wien2K. The resulting band structure is in excellent 
agreement to that from FPLO. Irrespective of the com- 
putational method, the LDA energy spectrum is metallic 
due to the underestimation of the correlation effects in 
this approximation. Experimentally, the green-colored 
BiCu2P06 is a magnetic insulator. The insulating be- 
havior is readily reproduced by the LSDA+C/ calculations 
(see Sec. HVBl . 

Eight Cu atoms in the crystallographic unit cell of 
BiCu2P06 give rise to eight 3d x 2_ y 2 bands (Fig. [3]). We 
first fit these bands with a tight-binding model and ex- 
tract the hopping parameters ti (Table H]) . The fitting 
procedure involves Wannier functions (WF) centered on 
Cu sites^ The application of the WF technique leads to 
a reliable fitting despite the slight overlap with the lower- 
lying bands. We are also able to resolve the couplings J 2 
and J 2 that correspond to the same Cu-Cu vector (0, 1, 0) 
but refer to different Cu sites in the structure (Fig. [I}. 
The hoppings are in agreement with the apparent fea- 
tures of the band structure. We find strong dispersion 
along the r — Y, T — Z, X — S, and U — R directions 
which correspond to the crystallographic fee-plane with 
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TABLE I. Leading hoppings of the tight-binding model and 
the resulting AFM exchange couplings. The exchange path- 
ways indicated in the first column are explicitly depicted in 
Figs. Q] and [4] The AFM part of the exchange integral is 
obtained by mapping the transfer integrals to an extended 
Hubbard model and eventually to a Heisenberg model using 
jAfm = ± t y UeS with UeS = 45 eV . 
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the couplings J\, J 2 , J 2 , J3, and J4. The dispersions 
along the other directions are less pronounced, indicative 
of a quasi-2D nature of this system. 

The hoppings are then introduced into a Hubbard 
model with the effective on-site Coulomb repulsion £/ e ff = 
4.5 eVi^Z r 29 ' 37 In the limit of strong correlations (ti <C 
U e ff) and in the half-filling regime, the low-lying excita- 
tions of the Hubbard model are described by a Heisen- 
berg Hamiltonian comprising AFM exchanges J AFM = 
Atf/U cS . The resulting J AFM values are listed in Ta- 
ble|H The maximum long-range hoppings ti beyond t\— 
amount to 30 meV, thus leading to J AFM < 10 K. Since 
the leading exchange couplings amount to 150 — 250 K, 
the minimal microscopic spin model can be restricted to 
five interactions: Ji, J2, J 2 , J3, and J4. 

A crucial fact to note at this juncture is the clear dif- 
ference in the strengths of J AFM and J.^™. Geometri- 
cally, the hopping paths for these exchanges are rather 
similar (Fig. Hjhand this structural feature led the au- 
thors of RefsTll5l and [H to assume J 2 — J 2 - In our 
analysis, we find that it is essential to treat these two ex- 
changes independently, otherwise the band splittings at 
the r point would not be reproduced correctly (i.e., one 
obtains four doubly-degenerate bands with J 2 = J 2 in- 
stead of the eight separate bands). Hence the frustrating 
next-nearest-neighbor exchanges "alternate" along the b 
axis (see Fig. [1]) with J 2 ~ O.5J2. A detailed analysis of 
this difference will be given in Sec. IIV CI 



B. LSDA+f/ 

The model approach allows to estimate all the ex- 
change couplings and to select the leading interactions 
for the minimum microscopic model. This is especially 
important for complex compounds with numerous and 
non-trivial superexchange pathways, like in BiC^POg. 
On the other hand, the model approach does not ac- 
count for FM contributions that are relevant for short- 
range interactions. 28 ' 38 To correct the leading couplings 
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for the FM contributions, we use the supercell approach. 
The total exchange integrals, consisting of the FM and 
AFM contributions, are listed in Table \H\ for the physi- 
cally reasonable range of the Ud values and for the two 
double-counting-correction (DCC) schemes. The latter 
is widely believed to be a minor feature of the LSDA+fJ 
method, but our recent studies evidenced a sizable influ- 
ence of the DCC on the exchange integrals in the case of 
short-range interactions £2r— 

The DCC is an essential part of the LSDA+J7 ap- 
proach, because a part of the on-site Coulomb repulsion 
energy is contained in LSDA and has to be subtracted 
from the total energy, after the explicit (mean-field) cor- 
rection for the on-site Coulomb repulsion is included. 
The two most common corrections are around-mean-ficld 
(AMF)ii and fully-localized-limit (FLL). 42 For spin-i 
magnetic insulators, the difference between AMF and 
FLL was commonly believed to be minor 42 By construc- 
tion, FLL looks more appropriate for the strongly local- 
ized regime U < U e s~& Yet, both AMF and FLL readily 
reproduce the insulating state of BiCu 2 P06. For exam- 
ple, we find the band gap E g ~ 2.4 eV and the magnetic 
moment of 0.81 fi B at U d = 6 eV in AMF. 45 FLL yields 
a somewhat lower gap E g ~ 1.6 eV at the same Ud value, 
but the gap is readily increased up to 2.1 eV at Ud = 8 eV. 
Experimental estimates of E g are presently lacking. How- 
ever, even the experimental input will hardly resolve the 
ambiguity, since the Ud value cannot be estimated pre- 
cisely. Then, the exchange couplings should be analyzed 
in more detail. 

AMF and FLL produce similar estimates for most of 
the couplings: Ji, J 2 , J 2 , and J 4 (see Table ITT]) . How- 
ever, the short-range interaction J3 is highly sensitive 
to the choice of the DCC. AMF suggests J3 to be a 
weak coupling (either FM or AFM, depending on the 
Ud value), while FLL ranks J3 as one of the leading 
AFM couplings, comparable to J% and J 2 . The FLL 
values essentially reproduce the previously published re- 
sults by Mentre et al^ that were also obtained within 
FLL but in a different band structure code. The model 
approach (TableU) evaluates J AFM , hence the FM contri- 
butions jf M — Ji — J AFM can be calculated. Following 
this procedure, we find a simple microscopic argument 
that supports the AMF results with weak J3. Both J\ 
and J3 arise from Cu-O-Cu superexchange with different 
angles at the oxygen atoms: 112.2° and 92.0°, respec- 
tively (see the top left panel of Fig.[4j. According to the 
Goodenough-Kanamori rules4£ the nearly 90° superex- 
change of J3 should yield the largest FM contribution. 
This conclusion conforms to the AMF results with J FM = 
-45 K and J FM = -135 K at U d = 6 eV. The FLL re- 
sults are opposite, jf M = -36 K and J FM = -16 K. As 
Ud is increased up to 8 eV, all the couplings are reduced, 
while the qualitative difference persists: |J FM | > |</ FM | 
in AMF, but I J 3 FM | < I jf M | in FLL. 

The above considerations suggest the exchange cou- 
plings from AMF as a more reliable estimate for 
BiCu 2 POg. For relevant examples from other compounds 



TABLE II. Total exchange couplings (in K) obtained from 
the LSDA+c^ calculations. The Ud value (in eV) denotes 
the Coulomb repulsion parameter of LSDA+c''. The last col- 
umn lists the double-counting correction scheme: around-the- 
mean-field (AMF) or fully-localized-limit (FLL). 
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with a simpler magnetic behavior, we refer the reader 
to Sec. IIV CI Additionally, we note that computational 
results for /3-Cu 2 V 2 7 (Ref. HJ and for several other 
Cu +2 -compounds^ also prefer AMF. Thus, we further 
rely on the AMF estimates and consider J3 as a weak 
coupling. The low value of J3 compared to J AFM re- 
duces the 2D Ji — J4 model, obtained from the model 
approach, to a quasi-lD model, depicted in the bottom 
part of Fig. Q] This model basically follows the earlier 
proposal by Mentre et alJ£ We find a two-leg spin lad- 
der with the leg coupling Ji, the rung coupling J4, and 
the next-nearest-neighbor frustrating couplings J 2 and J' 2 
along the legs. Yet, there are two important differences 
to be emphasized. First, the two next- nearest- neighbor 
couplings are inequivalent and fairly different. The J 2 
coupling connecting the Cu2 sites is twice as large as the 
coupling J 2 between the Cul sites (see Tables U and |H} . 
Second, we can safely establish the quasi-lD nature of the 
spin model, because the J3/J4 ratio is below 0.2 (com- 
pare to J3/ J4 = 0.55 — 0.65 in Ref. liH) . Both results are 
very important for understanding the material. 

The difference between J 2 and J 2 clearly alters the 
spin lattice. The pronounced one-dimensionality allows 
to simulate the behavior of the spin model on a quan- 
titative level, despite the presence of the strong frus- 
tration that narrows the range of applicable simulation 
techniques. Before turning to the experiments and sim- 
ulations (Sec. [V}, we will further discuss the non-trivial 
implementation of individual exchange couplings in the 
crystal structure of BiCu 2 P06 and provide further sup- 
port for the proposed spin model. 



C. Structural aspects of the magnetic exchange 

The interactions J\ and J3 run between corner-sharing 
and edge-sharing CUO4 plaquettes, respectively (top left 
panel of Fig. This geometry suggests Cu-O-Cu su- 
perexchange as the leading mechanism of the coupling 
and the angle at the oxygen atom as the key structural 
parameter determining the exchange integral. Follow- 
ing the Goodenough-Kanamori rules4^ we find that J3 
with the Cu-O-Cu angle of 92.0° is weakly AFM or 
even FM (see Table HI|) . The pathway of J\ reveals the 
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FIG. 4. (Color online) Parts of the crystal structure showing the details of individual superexchange pathways as well as the 
spin-ladder (top left panel) and the frustrated-spin-chain (right panel) features. The middle panel depicts the difference in 
the positions of the PO4 tetrahedra for the couplings J2 and J^. Curved arrows denote the rotations of the tetrahedra in the 
fictitious model structures (see text for details). The right panel shows the difference in the Ol-Ol distances for J2 and J2. 



sizably larger angle of 112.2° and, consequently, a siz- 
able AFM superexchange. A similar superexchange sce- 
nario is found in the mineral dioptase CueSigOis ■ 6H 2 
(green dioptase) 3 " 8 - and, presumably, in its anhydrous 
counterpart (black dioptase). The spin lattice of diop- 
tase comprises the AFM coupling J c between corner- 
sharing Cu04 plaquettes (the Cu-O-Cu angle amounts 
to 107.6° and 110.7° for green and black dioptase, re- 
spectively) and the FM coupling Jd between edge-sharing 
plaquettes (97.4° and 97.3°, respectively). More specif- 
ically, J c — 78 K and Jd = —37 K in green dioptase. 38 
The nature of the exchange couplings in the dioptase lat- 
tice is confirmed by the magnetic structure that was di- 
rectly investigated by neutron diffraction . 47 ' 48 Addition- 
ally, our recent computational study of green dioptase 
confirms the assignment of the exchange couplings and 
yields a consistent interpretation for all available exper- 
imental data^ 8 - The reference to the closely related su- 
perexchange scenario in dioptase should be taken as an 
additional argument for the weakness of J3 and the re- 
sulting quasi-lD character of BiCu 2 P06- 

In fact, one can find further examples supporting the 
pronounced difference between J\ and J3. Numerous 
cuprates with chains of edge-sharing plaquettes are ex- 
perimental realizations of frustrated spin chains with FM 
nearest-neighbor couplings. Such FM couplings arise 
from the Cu-O-Cu angle close to 90° and typically range 
from —100 K to —300 K for oxide compounds (e.g., 
Li 2 Cu0 2 , Li 2 CuZr0 4 )i^In BiCu 2 P0 6 , J 3 FM is smaller 
due to the folded arrangement of the plaquettes. Never- 
theless, the pronounced FM contribution reduces the to- 
tal exchange to a weak coupling, either FM or AFM, de- 
spite the sizable AFM contribution of J^ FM = 176 K (cf. 
Table[T]). The leg coupling J\ appears for the twisted con- 
figuration of corner-sharing plaquettes (see Fig. 2]) with 
the Cu-O-Cu angle of 112.2°. A similar configuration is 
found in AgCuVOi, where the angle amounts to 112.7°, 
and a pronounced AFM exchange coupling J ~ 300 K 
is founds Thus, our estimates of J\ and J3 are in line 
with the experience regarding other Cu compounds with 
firmly established microscopic models. 



All the above arguments support the quasi-lD model 
with weak J3. In the following, we use this model as a 
working hypothesis to interpret the magnetic behavior of 
BiCu 2 P06- The quasi-lD model captures the essential 
physics of the material, although certain features may 
require the extension of the model towards including J3 
or anisotropy effects (see Sec. I VII and IVIII| ) . 

Taking J3 as a weak interaction, we find J4 to be 
the leading coupling along the c direction. This cou- 
pling runs between the CUO4 plaquettes of neighbor- 
ing ribbons. The bonding between the ribbons arises 
from Bi cations (bottom left panel of Fig. @|, yet Bi 
does not give any sizable contribution to the states near 
the Fermi level. Therefore, we assign J4 to the Cu-O- 
O-Cu superexchange with the double 0-0 contact of 
2.75 A. Similar couplings between the disconnected cop- 
per plaquettes have been reported for (CuCl)LaNb 2 07 
and Bi 2 Cu04J 39 i 51 ' 52 Due to the large spatial separation 
of the Cu atoms (4.91 A), a sufficiently strong interac- 
tion arises for specific configurations of the ligand or- 
bitals only (see Ref. for an instructive example) . This 
explains the strong inter-ribbon coupling along the c di- 
rection, in contrast to a very weak coupling between the 
structural ribbons along a where the shortest Cu-Cu dis- 
tance is 4.85 A. 

Finally, we address the most puzzling feature of 
BiCu 2 P06, the next-nearest-neighbor couplings J 2 and 
J 2 . While the other couplings can be tentatively assigned 
after a careful analysis of the superexchange pathways, 
the sharp difference between J 2 and J 2 remains unex- 
pected. The Cu-Cu distances for the two couplings are 
the same and amount to 5.17 A, the lattice parameter 
along the b direction. On the other hand, J 2 and J 2 
correspond to different Cu positions and are inequiva- 
lent by symmetry. Band structure calculations within 
the model and LSDA+J7 approaches consistently suggest 
that jyj 2 ~ 0.5 (see Tables U and El) • 

The couplings J 2 and J 2 run between the copper pla- 
quettes, joined by another plaquette via Ol and by a 
PO4 tetrahedron via 02 (see the right panel of Fig. 0}. 
Thus, two different Cu-O-O-Cu channels are available. 



7 




& m ^ 

FIG. 5. (Color online) Wannier functions ("magnetic or- 
bitals" ) centered on Cu2 sites. Each orbital comprises the Cu 
3d x 2_ y 2 atomic orbital, large Ol and 02 <rp-contributions, 
and a smaller 03 <rp-contribution. 



Despite the very similar Cu-0 distances and Cu-O-0 
angles, there is a pronounced difference in the 01-01 
distances: 2.55 A for J 2 [the edge of the Cul plaque- 
tte] and 2.63 A for J' 2 [the edge of the Cu2 plaquette]. 
The shorter 01-01 distance should lead to the stronger 
coupling J2, in agreement with the computational result 
J 2 < J i- At first glance, the 02 channel looks completely 
identical, because the 02-02 distance is constrained by 
the edge of the PO4 tetrahedron (2.56 A). Nevertheless, 
this channel also contributes to the difference between J2 
and J' 2 . 

To get a deeper insight into the mechanism of the next- 
nearest-neighbor interactions, we inspect the Wannier 
functions for the Cul and Cu2 sites. Each WF comprises 
a Cu 3d x 2_ y 2 orbital along with the <r-type p-orbitals of 
the neighboring oxygens 01 and 02 (Fig. [5]). We also 
find small, but significant, u-contributions from second- 
neighbor oxygens 03 and 04 for the Cu2 and Cul WFs, 
respectively. These "tail" contributions arise from the 
specific orientation of the PO4 tetrahedra: one of the 
0-0 edges aligns along the Cu-02 bond, i.e., the Cu2- 
02-03 (<p) and Cul-02-04 (<p') angles approach 180°. 
Indeed, we find ip = 140.4° and ip' — 159.1° in agree- 
ment with the smaller 03 contribution of about 1.0 %, 
compared to 1.7 % for 04. 

Although the tail features of the WFs look tiny, they 
have a strong effect on the exchange couplings. To probe 
this, we constructed fictitious model structures by rotat- 
ing the PO4 tetrahedra around the 02-02 edge. Since 
the tetrahedra were kept rigid, only the ip and ip' an- 
gles were varied, while other geometrical parameters re- 
mained constant^ We found that the position of the 
tetrahedron leads to a dramatic change in the absolute 
values of J 2 and J 2 . As the ip angle is increased to- 
wards 180°, the 03 contribution gets larger, and J2 con- 
sequently decreases (Fig. [6]). The rotation of the tetra- 
hedra by 15° makes J2 and J' 2 equal, while the further 
rotation will switch the system to the J' 2 > J2 regime. 
The WF of Cul and the interaction J' 2 are less sensitive 

to the variation of the ip' angle within the studied angle 

54 

range. 

Our analysis shows that the structural features beyond 
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FIG. 6. (Color online) Exchange integrals J2 and J'2 and 
the contribution of the second- neighbor oxygens (03, 04) to 
the Wannier functions, depending on the position of the PO4 
tetrahedron. The dashed vertical line shows the angles in the 
BiCu2POe structure. 



the CUO4 plaquettes have a sizable effect on the exchange 
couplings in Cu +2 compounds. In BiCu2P06, the tails 
of the WFs on the second-neighbor oxygens have 90° ori- 
entation and should then reduce the AFM coupling (see 
Fig. EJ. This unexpected interference of the magnetic or- 
bitals on the second-neighbor oxygen site is one of the mi- 
croscopic reasons for the observed difference between J 2 
and J'2 ■ It is worth noting that the role of non-magnetic 
side groups was emphasized theoretically long ago^ but 
is often not taken into account adequately in a quanti- 
tative description. Here, we have shown that the oxygen 
orbitals play the key role, while the phosphorous atom 
simply "holds" the four oxygens of the tetrahedron to- 
gether. There is no appreciable phosphorous contribu- 
tion at the Fermi level, and its contribution to the WF's 
is also minor (below 0.1 %). This general mechanism, 
involving interacting oxygen atoms, has been recently 
found in vanadium phosphates^ and deserves further in- 
vestigation in the compounds comprising other transition 
metals. 



V. EXPERIMENTAL RESULTS 

A. Magnetic susceptibility 

The temperature dependence of the magnetic suscep- 
tibility is shown in Fig. [7] and resembles closely the data 
from Ref.[l5|. We find a broad maximum at T* ax ~ 62 K, 
indicative of the predominantly AFM low-dimensional 
and/or frustrated behavior. The sharp decrease in the 
susceptibility below Tj£ ax is a signature of the spin gap. 
In the low-temperature region, the 0.1 T data show a 
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FIG. 7. (Color online) Magnetic susceptibility of BiC^POs 
measured in the applied field /io-ff = 5 T and the fit of the ID 
spin model with g ~ 2.16, Ji ~ 140K, J2 = Ji, J'-z = § Ji, and 
J 4, — I Ji (simulation for a finite lattice with N = 20 sites). 
The inset shows the Curie- Weiss fit above 300 K. 

weak upturn below 5 K. This upturn is largely suppressed 
in the field of 5 T and can therefore be assigned to a para- 
magnetic contribution of defects/impurities. Above 10 K, 
the susceptibility is field-independent in the studied field 
range [IqH < 5 T. 

Above 200 K, the system approaches the Curie- Weiss 
regime. In order to improve previous studies; 15 ! 16 we 
measured the susceptibility at high temperatures up to 
700 K and fitted the data above 300 K with the expression 

where 6 is the Curie- Weiss temperature and C — 
N A (g^ B ) 2 S(S + l)/(3fcs) is the Curie constant. Our fit 
gives C = 0.447(1) emuK/mol Cu, and 9 = 181(1) K (see 
the inset of Fig. [7J- Fitting the data with an additional 
temperature- independent xo term leads to a small xo, 
therefore, we neglect this term in further analysis. We 
establish the predominant AFM nature of the exchange 
interactions with an energy scale of about 200 K. The C 
value corresponds to an effective moment of 1.89(1) \xb, 
slightly above the ideal spin-^ value of 1.73 /zg and rather 
typical for Cu+ 2 compounds^ Since T* ax /6> ~ 0.3, 
strong frustration should be expected. 

For further analysis, we fit the magnetic susceptibility 
using our microscopic spin model. Koteswararao et alJ^ 
have shown that the data do not conform to the model of 
isolated non-frustrated spin ladders. The introduction of 
interladder couplings does not significantly improve the 
description, 5 Therefore, realistic models with frustrating 
next-nearest-neighbor couplings have to be considered. 
Mentre et al*& used a frustrated J\ — Ji — J4 spin model 
and fitted the data with Ji ~ 140 K, J2 ~ 0.5Ji, and 
J 4 ~ 0.4 Ji (see also Ref. HH), but this model did not 
take into account the difference between J 2 and J' 2 . 

Here we employ the J\ — J2 — J' 2 — J a model (Fig. [TJ) to fit 



the experimental magnetic susceptibility. This ID frus- 
trated spin model can be treated by exact diagonaliza- 
tions for finite lattices or by renormalization-group tech- 
niques. The former turns out to be appropriate for the 
present problem due to the small finite-size effects and 
will be used here for the susceptibility fit. The unit cell 
comprises four inequivalent Cu 2+ ions, hence the num- 
ber of sites N in the finite cluster should be a multiple of 
four. To fit the experimental data, we first approximate 
our model by the following set of parameters Ji = J2, 
J 2 = \Ji, and J4 = |Ji, according to Table HU The 
simulations yield the reduced susceptibility x* which can 
be fitted to the experimentally observed x using 

X=^^X* (2) 

with only two variable parameters: g and J\. The sim- 
ulations for N = 16 and N — 20 sites provide almost 
identical susceptibility curves. Hence, finite-size effects 
are negligible and our simulations yield accurate results 
for the ID spin model under consideration. 

Our optimal fits yield Ji ~ 140 K and g ~ 2.16. 
The fitted g value is typical for Cu +2 compound a 29 i 58 
and also conforms to the effective magnetic moment of 
1.89 fiB which leads to g = 2.18. The absolute value 
of Ji is in remarkable agreement to the computational 
estimate of 100 - 150 K (cf. Table HJ). The fit follows 
the experimental data down to 100 K (see Fig. [7]). At 
lower temperatures, we find slight deviations from the 
experiment. For instance, the position of the suscepti- 
bility maximum T* ax is overestimated and the theoreti- 
cal curve lies slightly below the experimental data. This 
shows that our model overestimates the spin gap A. We 
shall return to this issue below in Sec. IVII 

We also tried to vary the ratios of exchange integrals 
and found several fits of similar quality. In particular, the 
parameter set from Ref.[l6l(Ji ~ 140 K, J 2 = J 2 ~ 0.5 Ji, 
and J4 ~ 0.4Ji) is also in agreement with the magnetic 
susceptibility data and yields a comparable g = 2.145. 
However, this parameter set does not account for the dif- 
ference between Ji and J 2 . Since, as shown above, this 
difference is evidenced by two different computational ap- 
proaches and has a clear structural origin, we regard the 
solution Ji = Ji, J' 2 = \Ji-, and J4 = § Ji as the micro- 
scopically justified parameter set for BiCu 2 P0 6 - 

B. High-field magnetization and the spin gap 

The low-energy physics of BiC^POg is characterized 
by the presence of a spin gap A. Previous estimates of 
A, based on the magnetic specific hea t 15 ' 17 and Knight 
shift consistently suggested A ~ 35 K. The INS data 
revealed a smaller gap of 2 meV (about 23 K)»i& The 
observed discrepancy calls for the application of further 
experimental methods, especially in light of the ambigu- 
ity of the specific heat and the Knight shift estimates, 
which arises from the fitting expressions that depend on 
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FIG. 8. (Color online) Magnetization curve of BiCu2POe 
measured in pulsed field at T = f.4 K. The magnetization 
values are given in arbitrary units (a.u.) The arrow shows 
the critical field H c where the spin gap is closed. The solid 
lines are guides for the eye. 



the character of the spin excitations and, in particular, 
on the dimensionality of the system. 

High-field magnetization data can provide a robust es- 
timate of the spin gap. The magnetization process of 
BiCu2PC>6 is presented in Fig. [5] At low fields, the mag- 
netization shows a weak linear increase with the field 
until fioH c ~ 22 T, where it bends upwards following a 
much steeper linear increase at higher fields^ The tran- 
sition at H c implies the closing of the spin gap and can be 
used for the numerical estimate of A. Similar to Ref. l3lL 
we take H c as the point of the maximum curvature. We 
find A = gfJ-B/J-oHc/kB — 32 K, in good agreement to the 
previous estimate A ~ 35 K obtained from the magnetic 
specific heat and the Knight shift data i 15 ' 17 ' 21 i 22 

The behavior of the magnetization for small fields 
needs to be discussed in more detail. In an ideal, SU(2) 
invariant and defect-free gapped system, the magnetiza- 
tion should be zero below H c (see also Fig. [T2l below) . 
Impurities give rise to a finite magnetization contribu- 
tion, but this should typically saturate around 5 T at 
the present low temperature of 1.4 K. Since the mea- 
sured magnetization keeps increasing up to 22 T, we con- 
clude that the weak linear field dependence for H < H c 
is due to the presence of weak anisotropic interactions 
in BiCi^POe- One such anisotropy, which is known^S to 
give rise to a linear magnetization response in the gapped 
regime of similar ladder systems, is the Dzyaloshinksy- 
Moriya (DM) anisotrop y 61 i 62 As explained in Ref. |6fJ, an 
isolated AFM dimer with a DM energy term of the form 
D • (Si x S2) admixes triplet excitations into the singlet 
ground state, and this gives rise to a uniform magnetiza- 
tion response of the form aDx(DxB) even far be- 
low the critical field. There is also a staggered response in 
first order in D of the form m s oc D x B which can be de- 
tected by a local probe, such as NMR experiments. Sim- 
ilar features arise in the spin-ladder Cu2(C5rIi2N 2 )2Cl4 



compoud. 60 Hence, it is reasonable to expect that the 
linear response observed for BiCu2P06 at H < H c stems 
from the presence of the DM anisotropy. 

For completeness, it is worth providing a brief dis- 
cussion on the main DM vectors, based on the crystal 
symmetry of BiCu 2 P0 6 (cf. Fig.Q}. First of all, a DM 
anisotropy on each rung is allowed by symmetry, since 
each rung comprises two inequivalent Cu sites and thus 
the inversion symmetry through the middle of each rung 
is lacking. The translational invariance along the b axis 
(with a period of two rungs) necessitates that the DM 
vectors are the same on every second rung. Furthermore, 
the fact that the ac plane is a reflection (i.e., crystallo- 
graphic mirror) planed confines the DM vectors to the b 
direction. There is finally a screw axis symmetry along 
b (translation along b by one rung, followed by a C2 ro- 
tation around the b axis) which connects the sites of two 
consecutive rungs. This last symmetry necessitates that 
the DM vectors on the two consecutive rungs differ in 
sign. The DM terms are also expected for other, inter- 
rung couplings. 

Finally, we would like to point out that the measured 
magnetization data right above H c do not show any 
square root singularity (cusp) as is typical for ID sys- 
tems with a quadratic branch of magnetic excitations 
above the ground state (see also Fig. IT21 below). This 
is probably related to the presence of the DM interac- 
tions mentioned above and the inter-ladder coupling J3, 
which arc both expected to smooth out the singularity. 

In Sec. I VIII below, we provide a more detailed theoret- 
ical picture for the magnetization process, but first it is 
essential to understand the nature of the lowest magnetic 
excitations in BiCu 2 P06- 



VI. LOW-ENERGY EXCITATIONS FROM 
EXACT DIAGONALIZATIONS 

We have performed an exact diagonalization study of 
the model Hamiltonian discussed above (but without DM 
terms) with parameters J± — J2 — 1, J 2 — 0.5, and 
J 4 = 0.75, using finite lattices of N = 12, 16, 20, 24, 28, 
32, and 36 sites with periodic boundary conditions along 
the legs (x-axis). The model is depicted in the upper 
panel of Fig. [5] Apart from translations along the legs 
by 2a, we also have two discrete spatial symmetries in 
this model. One is a reflection through any of the rungs 
(Vy) and the other is a 7r-rotation (£22) around the z- 
axis which is perpendicular to the plane of the ladder and 
passes through the center of a J\ — J4 rectangle. Instead 
of the latter, we can take the generator consisting of a 
translation by a, combined with a reflection along the 
axis crossing the middle of all rungs. 

A first strong insight into physics of this model comes 
from a simple examination of the ground state expecta- 
tion values of various local energy terms (sj • Sj). Owing 
to the spatial symmetries of the problem, there are four 
inequivalent bonds only. These are the bonds associated 
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FIG. 9. (Color online) Top panel: The actual structure (dis- 
regarding the buckling) of the present model, and in the lower 
panel its topologically equivalent version obtained by flipping 
the two sites of every second rung. 



with the four different exchange couplings J\ , J 2 , J' 2 , and 
J 4 in each unit cell. The corresponding ground state ex- 
pectation values, denoted as ex, &2, e' 2 , e 4 , are provided 
in Table IIIII together with the total ground state energy 
per site E/N = {2J x ex + J 2 e 2 + J 2 e' 2 + J 4 e 4 )/2. The 
latter shows only small finite-size variations for N > 16, 
which points to a very short correlation length^ More 
importantly, we observe a sizably large value for the spin- 
spin correlations on the rungs, e 4 ~ —0.47, which is more 
than twice the values on the remaining bonds. This re- 
sult tells us that the system is in the strong rung coupling 
regime, despite the fact that the leg couplings J\ and J 2 
are comparable to the rung coupling J4. 

In Fig. [TU1 we have superimposed the low-energy exci- 



TABLE III. The ground state expectation values of the four 
different bond strengths (sj ■ Sj) per unit cell and the total 
energy per site E/N in units of Ji. 



N 


ei 


e2 


e'2 


e 4 


E/N 


12 


-0.21568 


-0.17868 


-0.17162 


-0.42632 


-0.50780 


16 


-0.17466 


-0.19499 


-0.18411 


-0.47136 


-0.49494 


20 


-0.16633 


-0.21940 


-0.20928 


-0.46032 


-0.50096 


24 


-0.18409 


-0.18626 


-0.17460 


-0.47234 


-0.49800 


28 


-0.17210 


-0.20336 


-0.19237 


-0.47127 


-0.49860 


32 


-0.17683 


-0.19716 


-0.18571 


-0.47132 


-0.49858 



tations for each system size as a function of the allowed 
momentum quantum numbers so that we obtain a clear 
picture of the low-energy dispersion of the model. We 
observe that the lowest triplet (total spin 5=1) exci- 
tations (thick open symbols) form a well-defined (coher- 
ent) branch separated from the continuum by a finite gap 
for k > 0.47r/(2a). This branch has an incommensurate 
minimum at A: m i n ~ 0.87r/(2a) at A ED ~ 0.5Ji. In addi- 
tion to the lowest branch, we also find a second branch 
which is degenerate with the first at k = 7r/(2a) but this 
shifts quickly to higher energies into the continuum for 
k < 7r/(2a). 

Before we discuss the main implications of these results 
with regard to BiCu2P06, we would like to provide a ba- 
sic microscopic description of the excitation spectrum. 
To this end we perform a perturbative expansion around 
the limit of isolated rungs J\ — J 2 = J'2 = 0. We first 
introduce the singlet and triplet states of a single rung 
with sites 1 and 2 as |s) = (|t|) - |+t»/V2. |*i) = |tt), 
|t_i) = HI), and |*o) = (III) + lit)) A/2. The unper- 
turbed ground state is the product state of singlets on all 
rungs. Excitations arise by promoting one or more rungs 
into triplet states \t m ), with m = ±1,0. The inter-rung 
couplings have two effects. The first is that they renor- 
malize the ground state energy as well as the energies in 
the one-triplon sector. The second is that they induce 
a finite amplitude for nearest-neighbor and next-nearest- 
neighbor hoppings of triplons in the one-triplon sector. 
Including the amplitude from all different processes and 
exploiting the translational invariance by 2a, one finds 
two separate bands of one-triplon excitations due to the 
fact that we have two rungs per unit cell in the model. 
Their energies relative to the renormalized ground state 
energy are given by E K a ^(k) = A k ± \B k \, with 



Ak — J 4 



12J? + 3(J 2 + J'2) 2 - 4(J 2 - J 2 ) 



77 \2 



16J 4 

J2 + J'2 (J2 ~ J'2) ~ ^J\ 



8J 4 



cosfc 



M2 



(J2 + J'2 ) 

16J 4 



■ cos 2fc, 



jj fc = A(l + e -ifc)_ Jl ( J2 + J 2) 

2 V ' 8J 4 

x (l + e- lk + e Lk + e~ 2lk ) 

These second-order dispersions are shown in the lower 
panel of Fig. (TUJ Although its prediction for the spin gap 
is more than twice higher than the exact value (shown in 
the upper panel), the second-order perturbation theory 
captures well the position of the minimum and the overall 
shape of the dispersion. 

Next, we would like to comment that the degeneracy 
of the two branches at k = is not a generic feature 
of the exact dispersions but an accidental feature of the 
second-order expression for the given values of the ex- 
change integrals. In higher orders of perturbation theory 
or for slightly different parameter values, this degener- 
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FIG. 10. (Color online) Top: Superimposed low-energy dis- 
persions from exact diagonalizations on systems with N =12, 
16, 20, 24, 28, 32, and 36 sites and for the parameters 
Ji = J 2 = 1, J'i = Ji/2, and J 4 = 0.75 Ji. Empty (black) 
symbols denote the singlet S = states, thick open (blue) 
symbols denote the S = 1 states, and filled (red) symbols de- 
note the S — 2 states. The solid lines are polynomial fits to 
the visible parts of the lowest one-triplon excitation branches. 
Bottom: The two one-triplon energy branches predicted from 
second-order perturbation theory around the strong coupling 
limit (cf. text). We emphasize here that the degeneracy of 
the two bands at k = is an accidental feature of the second- 
order theory for the given values of the exchange parameters, 
while the degeneracy at k = 7r/(2a) is a generic feature re- 
lated to the fact that the model has a period a and not 2o 
along the legs of the ladder (cf. text). 



acy will be lifted. In contrast, the degeneracy of the two 
branches at k — n/(2a) is a generic feature and persists 
to all orders as seen in the exact spectra. The reason 
behind this is the presence of the discrete symmetry gen- 
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FIG. 11. (Color online) Top: Same as in Fig. HOI but in the 
symmetry setup of the lower panel of Fig. [9] The solid line is 
a polynomial fit to the visible part of the lowest one-triplon 
excitation branch. Bottom: The one-triplon energy dispersion 
predicted from second order perturbation theory around the 
strong coupling limit (cf. text). 



erator mentioned above (translation by a followed by a 
reflection through the middle of all rungs). To see this, 
we may start from the upper panel of Fig.|9]and exchange 
the two sites of every second rung without altering the 
topology of the model. This gives the equivalent model 
shown in the lower panel of Fig. [9] which has period a 
and not 2a. To elucidate this point, we may repeat the 
strong-coupling expansion in this alternative symmetry 
framework. To this end, we must take into account the 
extra negative signs that arise from the antisymmetry of 
the singlet rung wavefunction when nipping the two sites 
of every second rung. In terms of the new momenta, 
we now obtain a single one-triplon excitation band with 



energy dispersion 



E{k) = J 4 



12 J 2 + 3( J 2 + J' 2 ) 2 - 4( J 2 - J 2 ) 2 
16J 4 



c\ cos fc + C2 cos 2fc + C3 cos 3A; + C4 cos 4k (3) 



where ci = — J\ 
MJ2+4) 



Ji(.h+J' 2 ) „ _ J 2+ J ' 2 

4J 4 



C 2 2- 

(J 2 + J^) 2 



(J2-J^) 2 J? 
8J4 4J 4 ' 



c 3 = 'Vj ~ ) and C4 = — yJ2 1 ^ J j 2 ' ■ This dispersion is 
shown in the lower panel of Fig. (lTJ It is 



le lower panel of Fig. [TTJ It is clear that by 
folding this back into the Brillouin zone [— 7r/2a, it /2a] we 
shall obtain the two branches shown before in the lower 
panel of Fig. [TO] It is also evident in this representation 
that the incommensurate nature of the dispersion arises 
already in first order and is dictated by the frustrated 
couplings J2 and J' 2 which appear in the leading term in 
the above expression for c 2 . 

For completeness, we present the exact diagonalization 
results in the new symmetry setup in the upper panel of 
Fig. HU Again, the overall shape of the lowest dispersion 
and the position of the minimum are in agreement with 
the prediction of the strong coupling expansion shown in 
the lower panel. 

An interesting feature which becomes better visible in 
the representation of Fig. [IT] is the presence of a number 
of low- lying singlets for momenta close to k — ir/a. These 
excitations can be understood as a singlet bound state of 
two triplons. Such singlet excitations could be captured 
by optical experiments, such as phonon-assisted infrared 
absorption. Singlet bound states of two triplons have 
been observed using such techniques in cuprate ladders. 65 
In a broader context, the low- lying singlet at k = n/a can 
be considered as a singlet mode going soft at the transi- 
tion to a dimerized phase with dimers forming along the 
legs^ In the model considered, this scenario might occur 
as the rung coupling J4 is reduced further. 

Let us now discuss the implications of the above find- 
ings for BiCu 2 P0 6 . Taking J x ~ 140 K from the fit 
of the susceptibility we obtain for the spin gap A ED ~ 
0.5Ji ~ 70 K which is almost twice the value obtained 
from the high-field magnetization data, or the value re- 
ported by other groups i 15 ' 17 i 21 ' 22 Hence, we find that the 
present model of an isolated frustrated ladder overesti- 
mates the value of the spin gap in BiCu^POg, a fact 
that was already suggested from the behavior of the sus- 
ceptibility at low temperatures. One way to account for 
this discrepancy is to include a finite interladder coupling 
J3. Along the lines of the previous perturbative analy- 
sis, one finds that J3 gives rise to a first-order hopping 
of triplons along the y direction. As a result, the two 
bands attain a common extra dispersion term of the form 
— (J3/2) cos k y (with k y in units of tt divided by the inter- 
ladder distance). This shifts the minimum of the lowest 
band down by J3/2. Thus, to account for the 35 K spin 
gap one would need an interladder coupling of the order 
of J3 ~ 70K in this simple approximation. However, in 
the present regime we expect the perturbative calculation 
to be only qualitatively correct, so that a precise deter- 
mination of the interladder coupling either needs to come 
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FIG. 12. (Color online) Magnetization curve of BiCu2P06 as 
obtained from DMRG and ED. 



from more elaborate theoretical approaches (such as den- 
sity matrix renormalization group simulations of coupled 
ladders) or, ultimately, from inelastic neutron scattering 
experiments on single crystals. 



VII. MAGNETIZATION PROCESS FROM ED 
AND DMRG 

Here we revisit the magnetization process of 
BiCu2PC"6, in the light of the physical picture obtained 
above for the lowest magnetic excitations. To this end, 
we have employed Lanczos diagonalizations up to N — 
32 sites with periodic boundary conditions, as well as 
DMRG simulations with up to L = 128 rungs using open 
boundary conditions. Some representative magnetization 
curves are shown in Fig. 1121 The results from the two 
largest clusters treated by DMRG (L = 64, 128 rungs) 
converge to a rather smooth magnetization curve. They 
also give a critical field H c almost identical to the one ob- 
tained from ED for 32 sites, which further corroborates 
the value of the spin gap A ED ~ 0.5 Ji given above. 

To discuss the nature of the magnetization process 
in more detail, we distinguish three different regimes, 
namely the one at low magnetizations above H c , the one 
at high magnetizations as we approach the saturation 
field -ffsat, and the intermediate regime. The low mag- 
netization regime can be qualitatively understood on the 
basis of gradually filling the excitation band of Fig. QT] 
(bottom) with triplons as we ramp up the field above H c . 
One immediate consequence is the presence of a square- 
root singularity in the magnetization right above H c (cf. 
Fig. [T2|) which is due to the quadratic dispersion above 
the minimum. Another important ingredient in this con- 
sideration is the presence of two incommensurate minima 
(at k ~ ±0.47r/a) in the triplon dispersion which, given 
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FIG. 13. (Color online) One-magnon energy dispersions ob- 
tained analytically (see Eq. ©). 



the local hard-core constraint of the triplons, gives four 
Fermi points. Thus if the four-Fermi-point fix point is in- 
deed stable, the effective low-energy theory of BiCu2POg 
at low magnetizations is a two-component LL. 

In a similar way, the magnetization process close to 
the saturation field can be understood starting from the 
fully polarized state and gradually filling the one-magnon 
excitation bands by single spin flips. Using the setup of 
the lower panel of Fig. |9] and setting the energy of the 
fully polarized state to zero, one obtains two one-magnon 
bands which are given by the eigenvalues of the matrix 



H 



onc-magn 



v k u' k 



(4) 



with Uk = -(J4/2 + Ji + J 2 ) + J2 cos 2k, u' k = -( J 4 /2 + 
Ji + J 2 ) + J 2 cos2fc, and Vk = J4/2 + Jicosfc. The 
two one-magnon bands are shown in Fig. 1131 As ex- 
pected, we find that each band has two minima at in- 
commensurate wavevectors. In particular, the minima 
of the lowest band sit at k = ±0.4313l7r/a, which 
are close to the minimum fc-points of the triplon dis- 
persion of Fig. Qj] (bottom). The corresponding mini- 
mum energy E m i n = —3.5643 J\ gives the saturation field 
-Hsat = 3.5643 Ji/ (g/ioMs), hi agreement with the numer- 
ical results of Fig. [T5] By gradually filling the minimum 
of the lowest one-magnon branch, we describe the mag- 
netization process as we decrease the field below H sa ±. 
Similar to the low magnetization regime, the quadratic 
dispersion around the one-magnon minimum gives rise 
to a square-root singularity right below iJ sa t which can 
be seen in our numerical results of Fig. [T^l In addi- 
tion, the presence of two "incommensurate" minima (at 
±0.4313l7r/a) in the lowest one-magnon branch opens a 
possibility that the appropriate low-energy effective the- 
ory of BiCu2P06 at high fields is a two-component LL. 



It is presently unclear whether the possible two- 
component LL phases discussed at low and high fields 
form a single phase, or whether they are separated by 
one or more intervening phases at intermediate mag- 
netizations. Inspecting the numerical results displayed 
in Fig. 1121 a plateau might, for instance, occur at 
M = M sat /2. The phase immediately above the plateau 
also requires further investigation since there is a pos- 
sibility for a one-component LL phase before we reach 
the high-field two-component LL phase. Such a rich in- 
terplay between one- and two-component LL phases and 
plateaux is realized in the frustrated antiferromagnetic 
Ji — J2 Heisenberg chain model (see, e.g., Ref. |67| and 
references therein). Testing and confirming the scenario 
outlined here for the physics of BiCr^POg in high mag- 
netic fields requires a separate and more detailed investi- 
gation which is, however, beyond the scope of this article. 

Let us finally compare to the experimental magnetiza- 
tion data of Fig. [5J Given our earlier estimate of J\ ~ 
140 K from the susceptibility fit, we obtain H sat ~ 345 T, 
which is much larger than the range of fields accessible in 
our experiment (H max = 60 T). Hence, the highest mag- 
netization values reported in Fig. [5] correspond to less 
than 10% of M sat . In contrast to the above theoretical 
predictions, the measured magnetization does not show 
any square- root singularity right above H c . As we dis- 
cussed in Sec. IV Bl this gives evidence for inter-ladder 
coupling J3 and/or DM interactions which smooth out 
the singularity. 



VIII. DISCUSSION AND CONCLUSIONS 

Using DFT band structure calculations, we derived the 
minimum microscopic model of BiC^POv This model 
is based on a two-leg-ladder lattice and comprises four 
antiferromagnetic exchange couplings: J\ along the legs, 
J4 along the rungs, and the frustrating next-nearest- 
neighbor couplings J2 and J' 2 along the legs (Fig. [T]). Al- 
though such a model does not provide a complete and 
quantitative description of the compound, it is a rea- 
sonable compromise between the complexity of the sys- 
tem and the capabilities of present-day numerical simu- 
lation techniques for the evaluation of ground-state and 
finite-temperature properties of frustrated quantum spin 
systems. We showed that the ladder geometry leads to 
strong spin correlations on the rungs, despite the siz- 
able frustration and the weaker rung coupling. This 
feature might explain why the simple model of the un- 
frustrated spin ladder reproduces certain properties of 
BiCu 2 P06, especially the behavior upon the chemical 
substitution.— ~— On the other hand, the reduction to 
the simple ladder model cannot be justified microscop- 
ically, since the frustrating coupling J2 is of the same 
order as the leg and the rung couplings J\ and J4, re- 
spectively. In particular, the coupling Ji has an effect 
on the spin gap. The simple J\ — J4 two-leg ladder with 
J4 = I J\ shows a spin gap of about 0.3Ji^ while in our 
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model the gap amounts to 0.5Ji. Thus, the frustration 
enhances the gap in a spin ladder, similar to a conven- 
tional frustrated spin chain4 

We interpret BiCu2P06 as a system of two-leg ladders 
with frustrating couplings along the legs. The absolute 
values of individual exchange couplings leave an ambi- 
guity to describe the system as a frustrated spin ladder 
or as coupled frustrated spin chains. Indeed, the actual 
system shows features of both models. On the one hand, 
the strongest correlations are found on the rungs, as in 
ordinary ladders. On the other hand, the correlations 
along the legs are incommensurate and lead to the spin 
gap, being minimal at an incommensurate position in the 
Brillouin zone. 

BiCii2P06 is a peculiar spin-ladder system interesting 
for future investigation. One of the exciting branches 
could be high-field studies above H c . Recent experi- 
ments on (CsH^N^CuB^ evidenced the emergence of 
the LL physics in the high-field phase of the two-leg 
spin ladderi 11 ' 12 BiCu2P06 offers an opportunity to ex- 
plore similar effects in the presence of the frustration, 
where the incommensurate position of the gap might 
lead to a two-component LL or instabilities thereof at 
fields just above H c . Another advantage is the relative 
ease of the chemical substitution that has stimulated a 
range of experimental studies on Zn- and Ni-substituted 
samples^ - — Here, again, the incommensurate leg spin- 
spin correlations could influence the effective interaction 
mediated between the impurity-induced localized spins, 
and thus lead to hitherto unobserved frustration effects. 

While working on the minimal microscopic model, one 
also has to understand its limitations. The main and 
most severe limitation is the reduction to a purely ID 
regime by neglecting J3. In fact, our band structure cal- 
culations suggest J3/J4 < 0.2, i.e., \J$\ < 25 K. If we 
adjust J3 to account for the actual spin gap A ~ 32 K 
~ 0.2 Ji, a larger value is obtained (see Sec. IVII) . Addi- 
tionally, the shape of the magnetization curve with the 
linear increase right above H c (Sec. IV B[) may exclude 
a purely ID scenario and point to sizable inter-ladder 
couplings. Considering all these arguments, we conclude 
that the inter-ladder coupling J 3 is likely relevant for the 
full picture, but its accurate estimate remains a chal- 
lenging task. Band structure calculations equally allow 
for FM or AFM J 3 (Table HE}- Experimental estimates 
would require theoretical information on a complex 2D 
Ji — J2 — J2 — J3 — J a frustrated spin system with long- 
range couplings Ji and J^. Such a system is basically 
beyond the capabilities of present-day numerical meth- 
ods. Therefore, the most reasonable approach could be 
analytical perturbation treatment, based on the accu- 
rate results for the ID model. We believe that this ap- 
proach will help to clarify the complex magnetic behav- 
ior of BiCr^POg and to improve the theoretical estimate 
of the spin gap with respect to the experimental value 
A ~ 32 K. 



The second limitation of our model is the lack of 
anisotropy effects. In particular, the DM interactions 
scale with J and can be sizable due to the strong isotropic 
exchange of 100 — 150 K. The DM couplings are al- 
lowed for all the bonds of the spin lattice with few re- 
strictions on the arrangement of the D vectors with re- 
spect to the crystal axes (see also Sec. IV Bj) . The com- 
prehensive investigation of the anisotropy effects would 
require electron spin resonance measurements on single 
crystals along with sophisticated band structure calcula- 
tions. Presently, we note that the increase in the mag- 
netization below H c (Fig. [5]) is a possible signature of 
the DM couplings. The non-zero Knight shift at low 
temperature s 21 ' 23 may have the same origin. 

In summary, our study provides a comprehensive de- 
scription of isotropic exchange couplings in the spin-i 
quantum magnet BiC^POg. We interpret this com- 
pound as a two-leg spin ladder with frustrating next- 
nearest-neighbor couplings along the legs. The leg cou- 
pling (Ji), the rung coupling (J4), and one of the next- 
nearest-neighbor couplings {J 2) amount to 120 — 150 K, 
while the other next-nearest-neighbor coupling J' 2 is half 
of J 2 due to the subtle structural differences between 
the respective superexchange pathways. The complex 
crystal structure of the compound leads to a non-trivial 
implementation of the spin ladder with two legs resid- 
ing on different structural ribbons. The proposed spin 
model is a derivative of the simple two-leg spin ladder 
and shows leading spin correlations on the rungs. Frus- 
trating couplings increase the spin gap and induce the 
incommensurate minimum of the triplon dispersion as 
well as an exotic behavior in high magnetic fields. The 
effects beyond our spin model include the inter-ladder 
coupling and the anisotropy. Experimental data show 
possible signatures of these effects and call for further 
investigation of BiCu2POg by means of inelastic neutron 
scattering and electron spin resonance measurements on 
single crystals. 
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